Load the necessary libraries
library(rstanarm) # for fitting models in STAN
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(tidyverse) # for data wrangling etc
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
theme_set(theme_grey()) # put the default ggplot theme back
Here is an example from Fowler, Cohen, and Jarvis (1998). An agriculturalist was interested in the effects of fertilizer load on the yield of grass. Grass seed was sown uniformly over an area and different quantities of commercial fertilizer were applied to each of ten 1 m2 randomly located plots. Two months later the grass from each plot was harvested, dried and weighed. The data are in the file fertilizer.csv in the data folder.
| FERTILIZER | YIELD |
|---|---|
| 25 | 84 |
| 50 | 80 |
| 75 | 90 |
| 100 | 154 |
| 125 | 148 |
| ... | ... |
| FERTILIZER: | Mass of fertilizer (g.m-2) - Predictor variable |
| YIELD: | Yield of grass (g.m-2) - Response variable |
The aim of the analysis is to investigate the relationship between fertilizer concentration and grass yield.
We will start off by reading in the Fertilizer data. There are many functions in R that can read in a CSV file. We will use a the read_csv() function as it is part of the tidyverse ecosystem.
fert <- read_csv("../public/data/fertilizer.csv", trim_ws = TRUE)
glimpse(fert)
## Rows: 10
## Columns: 2
## $ FERTILIZER <dbl> 25, 50, 75, 100, 125, 150, 175, 200, 225, 250
## $ YIELD <dbl> 84, 80, 90, 154, 148, 169, 206, 244, 212, 248
The individual responses (\(y_i\), observed yields) are each expected to have been independently drawn from normal (Gaussian) distributions (\(\mathcal{N}\)). These distributions represent all the possible yields we could have obtained at the specific (\(i^th\)) level of Fertilizer. Hence the \(i^th\) yield observation is expected to have been drawn from a normal distribution with a mean of \(\mu_i\).
Although each distribution is expected to come from populations that differ in their means, we assume that all of these distributions have the same variance (\(\sigma^2\)).
We need to supply priors for each of the parameters to be estimated (\(\beta_0\), \(\beta_1\) and \(\sigma\)). Whilst we want these priors to be sufficiently vague as to not influence the outcomes of the analysis (and thus be equivalent to the frequentist analysis), we do not want the priors to be so vague (wide) that they permit the MCMC sampler to drift off into parameter space that is both illogical as well as numerically awkward.
As a starting point, lets assign the following priors:
Note, when fitting models through either rstanarm or brms, the priors assume that the predictor(s) have been centred and are to be applied on the link scale. In this case the link scale is an identity.
Model formula: \[ \begin{align} y_i &\sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i &= \beta_0 + \beta_1 x_i\\ \beta_0 &\sim{} \mathcal{N}(0,100)\\ \beta_1 &\sim{} \mathcal{N}(0,10)\\ \sigma &\sim{} \mathcal{cauchy}(0,5)\\ OR\\ \sigma &\sim{} \mathcal{Exp}(1)\\ OR\\ \sigma &\sim{} \mathcal{gamma}(2,1)\\ \end{align} \]
Note, for routines that take more than a couple of seconds to perform (such as most Bayesian models), it is a good idea to cache the Rmarkdown chunks in which the routine is performed. That way, the routine will only be run the first time and any objects generated will be stored for future use. Thereafter, provided the code has not changed, the routine will not be re-run. Rather, knitr will just retrieve the cached objects and continue on.
For the purpose of comparing the frequentist and Bayesian outcomes, it might be useful to start by fitting the frequentist simple linear model. We will also fit this model with the predictor (fertilizer concentration) centred.
summary(lm(YIELD ~ FERTILIZER, data = fert))
##
## Call:
## lm(formula = YIELD ~ FERTILIZER, data = fert)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.79 -11.07 -5.00 12.00 29.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.93333 12.97904 4.001 0.00394 **
## FERTILIZER 0.81139 0.08367 9.697 1.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19 on 8 degrees of freedom
## Multiple R-squared: 0.9216, Adjusted R-squared: 0.9118
## F-statistic: 94.04 on 1 and 8 DF, p-value: 1.067e-05
summary(lm(YIELD ~ scale(FERTILIZER, scale = FALSE), data = fert))
##
## Call:
## lm(formula = YIELD ~ scale(FERTILIZER, scale = FALSE), data = fert)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.79 -11.07 -5.00 12.00 29.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 163.50000 6.00813 27.213 3.58e-09 ***
## scale(FERTILIZER, scale = FALSE) 0.81139 0.08367 9.697 1.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19 on 8 degrees of freedom
## Multiple R-squared: 0.9216, Adjusted R-squared: 0.9118
## F-statistic: 94.04 on 1 and 8 DF, p-value: 1.067e-05
In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.
fert.rstanarm <- stan_glm(YIELD ~ FERTILIZER,
data = fert,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
In the above:
Having allowed rstanarm to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the rstanarm development team make it very clear that there is no guarantee that the default priors will remain the same into the future.
prior_summary(fert.rstanarm)
## Priors for model 'fert.rstanarm'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 164, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 164, scale = 160)
##
## Coefficients
## Specified prior:
## ~ normal(location = 0, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 0, scale = 2.1)
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.016)
## ------
## See help('prior_summary.stanreg') for more details
This tells us:
mean(fert$YIELD)
## [1] 163.5
sd(fert$YIELD)
## [1] 63.97439
2.5 * sd(fert$YIELD)
## [1] 159.936
2.5 * sd(fert$YIELD) / sd(fert$FERTILIZER)
## [1] 2.113004
rstanarm, this is a exponential with a rate of 1 which is then adjusted by division with the standard deviation of the response.1 / sd(fert$FERTILIZER)
## [1] 0.01321157
One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging predictions would ensure that the priors are unlikely to influence the actual predictions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of parameter space that are not supported by the actual data, the sampler can become unstable and have difficulty.
We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refitting the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.
fert.rstanarm1 <- update(fert.rstanarm, prior_PD = TRUE)
ggpredict(fert.rstanarm1) %>% plot(add.data = TRUE)
## $FERTILIZER
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
I will also overlay the raw data for comparison.
fert.rstanarm2 <- stan_glm(YIELD ~ FERTILIZER,
data = fert,
prior_intercept = normal(164, 10, autoscale = FALSE),
prior = normal(0, 1, autoscale = FALSE),
prior_aux = cauchy(0, 2, autoscale = FALSE),
prior_PD = TRUE,
iter = 5000, warmup = 1000,
chains = 3, thin = 5, refresh = 0
)
ggpredict(fert.rstanarm2) %>%
plot(add.data = TRUE)
## $FERTILIZER
Now lets refit, conditioning on the data.
fert.rstanarm3 <- update(fert.rstanarm2, prior_PD = FALSE)
posterior_vs_prior(fert.rstanarm3,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y")
)
Conclusions:
ggpredict(fert.rstanarm3) %>% plot(add.data = TRUE)
## $FERTILIZER
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
fert.brm <- brm(bf(YIELD ~ FERTILIZER),
data = fert,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
In the above:
brms can define models that are not possible by most other routines. To facilitate this enhanced functionality, we usually define a brms formula within its own bf() function along with the family (in this case, it is Gaussian, which is the default and therefore can be omitted.)The returned object (a list) contains a range of objects. The output from rstan is contained in the fit list item.
fert.brm %>% names()
## [1] "formula" "data" "prior" "data2" "stanvars" "model"
## [7] "ranef" "save_pars" "algorithm" "backend" "threads" "fit"
## [13] "criteria" "file" "version" "family" "autocor" "cov_ranef"
## [19] "stan_funs" "data.name"
Having allowed brms to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the brms development team make it very clear that there is no guarantee that the default priors will remain the same into the future.
prior_summary(fert.brm)
## prior class coef group resp dpar nlpar bound source
## (flat) b default
## (flat) b FERTILIZER (vectorized)
## student_t(3, 161.5, 90.4) Intercept default
## student_t(3, 0, 90.4) sigma default
This tells us:
median(fert$YIELD)
## [1] 161.5
mad(fert$YIELD)
## [1] 90.4386
mad is the median absolute deviation. mad is to var as median is to mean.
standist::visualize("student_t(3,161.5,90.4)", xlim = c(-100, 1000))
for the beta coefficients (in this case, just the slope), the default prior is a improper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.
the prior on sigma is also a student t (although only the positive half is used in the stan code), with mean of 0 and standard deviation of 90.4.
standist::visualize("student_t(3,0,90.4)", xlim = c(-10, 500))
One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging predictions would ensure that the priors are unlikely to influence the actual predictions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of parameter space that are not supported by the actual data, the sampler can become unstable and have difficulty.
In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior = 'only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.
Lets try a Gaussian (normal) distribution for the slope. We will start with what is likely to be a very wide distribution (wide with respect to what we know about the likely slope).
standist::visualize("normal(0, 10)", xlim = c(-50, 50))
standist::visualize("normal(0, 0.1)", xlim = c(-1, 1))
fert.brm1 <- brm(bf(YIELD ~ FERTILIZER),
data = fert,
prior = prior(normal(0, 10), class = "b"),
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
fert.brm1 %>%
ggpredict() %>%
plot(add.data = TRUE)
## $FERTILIZER
fert.brm1 %>%
ggemmeans(~FERTILIZER) %>%
plot(add.data = TRUE)
fert.brm1 %>%
conditional_effects() %>%
plot(points = TRUE)
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (which I have mainly plucked out of thin air):
standist::visualize("normal(164, 10)", xlim = c(-100, 300))
standist::visualize("normal(0, 1)", xlim = c(-10, 10))
standist::visualize("cauchy(0, 2)", xlim = c(-10, 10))
standist::visualize("gamma(2, 1)", xlim = c(0, 10))
standist::visualize(
"student_t(3, 0, 90.4)",
"cauchy(0, 2)",
"gamma(2, 1)",
"gamma(2, 0.5)",
"gamma(2, 0.2)",
xlim = c(0, 20)
)
I will also overlay the raw data for comparison.
priors <- prior(normal(164, 10), class = "Intercept") +
prior(normal(0, 1), class = "b") +
prior(gamma(2, 1), class = "sigma")
fert.brm2 <- brm(bf(YIELD ~ FERTILIZER),
data = fert,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)
fert.brm2 %>%
ggemmeans(~FERTILIZER) %>%
plot(add.data = TRUE)
# OR
fert.brm2 %>%
conditional_effects() %>%
plot(points = TRUE)
fert.brm3 <- update(fert.brm2, sample_prior = "yes", refresh = 0)
When we have indicated that the posterior should be informed by both the prior and the posterior, both prior (governed by priors alone) and posterior (governed by both priors and data/likelihood) draws are returned. These can be compared by exploring the probabilities associated with specific hypotheses - the most obvious of which is that of no effect (that the parameter = 0).
When doing so, ideally the posterior should be distinct from the prior. If this is not the case, it might suggest that the posteriors are being too strongly driven by the priors.
fert.brm3 %>% get_variables()
## [1] "b_Intercept" "b_FERTILIZER" "sigma" "prior_Intercept"
## [5] "prior_b" "prior_sigma" "lp__" "accept_stat__"
## [9] "stepsize__" "treedepth__" "n_leapfrog__" "divergent__"
## [13] "energy__"
## fert.brm3 %>% hypothesis('Intercept=0', class='b') %>% plot
fert.brm3 %>%
hypothesis("FERTILIZER = 0") %>%
plot()
## fert.brm3 %>% hypothesis('scaleFERTILIZERscaleEQFALSE=0') %>% plot
fert.brm3 %>%
hypothesis("sigma = 0", class = "") %>%
plot()
Unfortunately, it is not possible to do this comparison sensibly for the intercept. The reason for this is that the prior for intercept was applied to an intercept that is associated with centred continuous predictors (predictors are centred behind the scenes). Since we did not centre the predictor, the intercept returned is as if uncentred. Hence, the prior and posterior are on different scales.
fert.brm3 %>% get_variables()
## [1] "b_Intercept" "b_FERTILIZER" "sigma" "prior_Intercept"
## [5] "prior_b" "prior_sigma" "lp__" "accept_stat__"
## [9] "stepsize__" "treedepth__" "n_leapfrog__" "divergent__"
## [13] "energy__"
fert.brm3 %>%
posterior_samples() %>%
select(-`lp__`) %>%
pivot_longer(cols = everything(), names_to = "key", values_to = "value") %>%
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "b"),
Class = ifelse(str_detect(key, "Intercept"), "Intercept",
ifelse(str_detect(key, "sigma"), "Sigma", "b")
)
) %>%
ggplot(aes(x = Type, y = value)) +
stat_pointinterval() +
facet_wrap(~Class, scales = "free")
fert.brm3 %>% standata()
## $N
## [1] 10
##
## $Y
## [1] 84 80 90 154 148 169 206 244 212 248
##
## $K
## [1] 2
##
## $X
## Intercept FERTILIZER
## 1 1 25
## 2 1 50
## 3 1 75
## 4 1 100
## 5 1 125
## 6 1 150
## 7 1 175
## 8 1 200
## 9 1 225
## 10 1 250
## attr(,"assign")
## [1] 0 1
##
## $prior_only
## [1] 0
##
## attr(,"class")
## [1] "standata"
fert.brm3 %>% stancode()
## // generated with brms 2.14.4
## functions {
## }
## data {
## int<lower=1> N; // total number of observations
## vector[N] Y; // response variable
## int<lower=1> K; // number of population-level effects
## matrix[N, K] X; // population-level design matrix
## int prior_only; // should the likelihood be ignored?
## }
## transformed data {
## int Kc = K - 1;
## matrix[N, Kc] Xc; // centered version of X without an intercept
## vector[Kc] means_X; // column means of X before centering
## for (i in 2:K) {
## means_X[i - 1] = mean(X[, i]);
## Xc[, i - 1] = X[, i] - means_X[i - 1];
## }
## }
## parameters {
## vector[Kc] b; // population-level effects
## real Intercept; // temporary intercept for centered predictors
## real<lower=0> sigma; // residual SD
## }
## transformed parameters {
## }
## model {
## // likelihood including all constants
## if (!prior_only) {
## target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
## }
## // priors including all constants
## target += normal_lpdf(b | 0, 1);
## target += normal_lpdf(Intercept | 164, 10);
## target += gamma_lpdf(sigma | 2, 1);
## }
## generated quantities {
## // actual population-level intercept
## real b_Intercept = Intercept - dot_product(means_X, b);
## // additionally draw samples from priors
## real prior_b = normal_rng(0,1);
## real prior_Intercept = normal_rng(164,10);
## real prior_sigma = gamma_rng(2,1);
## // use rejection sampling for truncated priors
## while (prior_sigma < 0) {
## prior_sigma = gamma_rng(2,1);
## }
## }
MCMC sampling behaviour
available_mcmc()
| Package | Description | function | rstanarm | brms |
|---|---|---|---|---|
| bayesplot | Traceplot | mcmc_trace |
plot(mod, plotfun='trace') |
mcmc_plot(mod, type='trace') |
| Density plot | mcmc_dens |
plot(mod, plotfun='dens') |
mcmc_plot(mod, type='dens') |
|
| Density & Trace | mcmc_combo |
plot(mod, plotfun='combo') |
mcmc_plot(mod, type='combo') |
|
| ACF | mcmc_acf_bar |
plot(mod, plotfun='acf_bar') |
mcmc_plot(mod, type='acf_bar') |
|
| Rhat hist | mcmc_rhat_hist |
plot(mod, plotfun='rhat_hist') |
mcmc_plot(mod, type='rhat_hist') |
|
| No. Effective | mcmc_neff_hist |
plot(mod, plotfun='neff_hist') |
mcmc_plot(mod, type='neff_hist') |
|
| rstan | Traceplot | stan_trace |
stan_trace(mod) |
stan_trace(mod) |
| ACF | stan_ac |
stan_ac(mod) |
stan_ac(mod) |
|
| Rhat | stan_rhat |
stan_rhat(mod) |
stan_rhat(mod) |
|
| No. Effective | stan_ess |
stan_ess(mod) |
stan_ess(mod) |
|
| Density plot | stan_dens |
stan_dens(mod) |
stan_dens(mod) |
|
| ggmcmc | Traceplot | ggs_traceplot |
ggs_traceplot(ggs(mod)) |
ggs_traceplot(ggs(mod)) |
| ACF | ggs_autocorrelation |
ggs_autocorrelation(ggs(mod)) |
ggs_autocorrelation(ggs(mod)) |
|
| Rhat | ggs_Rhat |
ggs_Rhat(ggs(mod)) |
ggs_Rhat(ggs(mod)) |
|
| No. Effective | ggs_effective |
ggs_effective(ggs(mod)) |
ggs_effective(ggs(mod)) |
|
| Cross correlation | ggs_crosscorrelation |
ggs_crosscorrelation(ggs(mod)) |
ggs_crosscorrelation(ggs(mod)) |
|
| Scale reduction | ggs_grb |
ggs_grb(ggs(mod)) |
ggs_grb(ggs(mod)) |
|
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
plot(fert.rstanarm3, plotfun = "mcmc_trace")
The chains appear well mixed and very similar
plot(fert.rstanarm3, "acf_bar")
There is no evidence of autocorrelation in the MCMC samples
plot(fert.rstanarm3, "rhat_hist")
All Rhat values are below 1.05, suggesting the chains have converged.
neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
plot(fert.rstanarm3, "neff_hist")
Ratios all very high.
plot(fert.rstanarm3, "combo")
plot(fert.rstanarm3, "violin")
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
stan_trace(fert.rstanarm3)
The chains appear well mixed and very similar
stan_ac(fert.rstanarm3)
There is no evidence of autocorrelation in the MCMC samples
stan_rhat(fert.rstanarm3)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
stan_ess(fert.rstanarm3)
Ratios all very high.
stan_dens(fert.rstanarm3, separate_chains = TRUE)
The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
fert.ggs <- ggs(fert.rstanarm3)
ggs_traceplot(fert.ggs)
The chains appear well mixed and very similar
ggs_autocorrelation(fert.ggs)
There is no evidence of autocorrelation in the MCMC samples
ggs_Rhat(fert.ggs)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
ggs_effective(fert.ggs)
Ratios all very high.
ggs_crosscorrelation(fert.ggs)
ggs_grb(fert.ggs)
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
fert.brm3 %>% mcmc_plot(type = "combo")
fert.brm3 %>% mcmc_plot(type = "trace")
fert.brm3 %>% mcmc_plot(type = "dens_overlay")
The chains appear well mixed and very similar
fert.brm3 %>% mcmc_plot(type = "acf_bar")
There is no evidence of autocorrelation in the MCMC samples
fert.brm3 %>% mcmc_plot(type = "rhat_hist")
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
fert.brm3 %>% mcmc_plot(type = "neff_hist")
Ratios all very high.
fert.brm3 %>% mcmc_plot(type = "combo")
fert.brm3 %>% mcmc_plot(type = "violin")
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
fert.brm3$fit %>% stan_trace()
fert.brm3$fit %>% stan_trace(inc_warmup = TRUE)
The chains appear well mixed and very similar
fert.brm3$fit %>% stan_ac()
There is no evidence of autocorrelation in the MCMC samples
fert.brm3$fit %>% stan_rhat()
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
fert.brm3$fit %>% stan_ess()
Ratios all very high.
fert.brm3$fit %>% stan_dens(separate_chains = TRUE)
The ggmcmc package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
fert.ggs <- fert.brm3$fit %>% ggs(inc_warmup = TRUE, burnin = TRUE) # does not seem to ignore burnin?
fert.ggs %>% ggs_traceplot()
fert.ggs <- fert.brm3$fit %>% ggs(inc_warmup = FALSE, burnin = TRUE) # does not seem to ignore burnin?
fert.ggs %>% ggs_traceplot()
The chains appear well mixed and very similar
fert.ggs %>% ggs_autocorrelation()
There is no evidence of autocorrelation in the MCMC samples
fert.ggs %>% ggs_Rhat()
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
fert.ggs %>% ggs_effective()
Ratios all very high.
fert.ggs %>% ggs_crosscorrelation()
fert.ggs %>% ggs_grb()
Posterior probability checks
available_ppc()
| Package | Description | function | rstanarm | brms |
|---|---|---|---|---|
| bayesplot | Density overlay | ppc_dens_overlay |
pp_check(mod, plotfun='dens_overlay') |
pp_check(mod, type='dens_overlay') |
| Obs vs Pred error | ppc_error_scatter_avg |
pp_check(mod, plotfun='error_scatter_avg') |
pp_check(mod, type='error_scatter_avg') |
|
| Pred error vs x | ppc_error_scatter_avg_vs_x |
pp_check(mod, x=, plotfun='error_scatter_avg_vs_x') |
pp_check(mod, x=, type='error_scatter_avg_vs_x') |
|
| Preds vs x | ppc_intervals |
pp_check(mod, x=, plotfun='intervals') |
pp_check(mod, x=, type='intervals') |
|
| Partial plot | ppc_ribbon |
pp_check(mod, x=, plotfun='ribbon') |
pp_check(mod, x=, type='ribbon') |
|
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
pp_check(fert.rstanarm3, plotfun = "dens_overlay")
The model draws appear to be consistent with the observed data.
pp_check(fert.rstanarm3, plotfun = "error_scatter_avg")
There is no obvious pattern in the residuals.
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "error_scatter_avg_vs_x")
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "intervals")
pp_check(fert.rstanarm3, x = fert$FERTILIZER, plotfun = "ribbon")
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
# library(shinystan)
# launch_shinystan(fert.rstanarm3)
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
preds <- posterior_predict(fert.rstanarm3, nsamples = 250, summary = FALSE)
fert.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = fert$YIELD,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = "gaussian"
)
plot(fert.resids)
Conclusions:
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
fert.brm3 %>% pp_check(type = "dens_overlay", nsamples = 100)
The model draws appear to be consistent with the observed data.
fert.brm3 %>% pp_check(type = "error_scatter_avg")
There is no obvious pattern in the residuals.
fert.brm3 %>% pp_check(x = "FERTILIZER", type = "error_scatter_avg_vs_x")
fert.brm3 %>% pp_check(x = "FERTILIZER", type = "intervals")
pp_check(fert.brm3, x = "FERTILIZER", type = "ribbon")
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
# library(shinystan)
# launch_shinystan(fert.brm3)
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
preds <- fert.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
fert.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = fert$YIELD,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE
)
fert.resids %>% plot()
The integerResponse argument indicates that noise will be added to the residuals such that they will become integers. This is important in for Poisson, Negative Binomial and Binomial models.
Conclusions:
fert.rstanarm3 %>%
ggpredict() %>%
plot(add.data = TRUE)
## $FERTILIZER
fert.rstanarm3 %>%
ggemmeans(~FERTILIZER) %>%
plot(add.data = TRUE)
fert.rstanarm3 %>%
fitted_draws(newdata = fert) %>%
median_hdci() %>%
ggplot(aes(x = FERTILIZER, y = .value)) +
geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
geom_line()
fert.brm3 %>%
ggemmeans(~FERTILIZER) %>%
plot(add.data = TRUE)
fert.brm3 %>% conditional_effects()
fert.brm3 %>%
conditional_effects() %>%
plot(points = TRUE)
fert.brm3 %>%
conditional_effects(spaghetti = TRUE, nsamples = 200) %>%
plot(points = TRUE)
fert.brm3 %>%
ggpredict() %>%
plot(add.data = TRUE)
## $FERTILIZER
fert.brm3 %>%
ggemmeans(~FERTILIZER) %>%
plot(add.data = TRUE)
fert.brm3 %>%
fitted_draws(newdata = fert) %>%
median_hdci() %>%
ggplot(aes(x = FERTILIZER, y = .value)) +
geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
geom_line()
| Package | Function | Description |
|---|---|---|
as.matrix() |
Returns \(n\times p\) matrix | |
tidybayes |
fitted_draws(newdata) |
Returns \(n\times p\) tibble with additional info about the chains, iteration and draw. Predictions of newdata |
tidybayes |
linpred_draws() | Returns $n\times p$ tibble with additional info about the chains, interaction and draw. Predictions on link scale | |tidybayes|tidy_draws()| Returns $n\times p$ tibble with addition info about the chain, iteration and draw | |tidybayes|spread_draws()| Returns $n\times r$ tibble (where $r$ is the number of requested parameters) with additional info about chain, interaction and draw | |tidybayes|gather_draws()| Returns a gatheredspread_drawstibble with additional info about chain, interaction and draw | |brms|posterior_samples()` |
Returns \(n\times p\) data.frame |
where \(n\) is the number of MCMC samples and \(p\) is the number of parameters to estimate. For the tidybayes versions in the table above, the function expects a model to be the first parameter (and a dataframe to be the second). There are also add_ versions which expect a dataframe to be the first argument and the model to be the second. These alternatives facilitate pipings with different starting objects.
| Function | Description |
|---|---|
median_qi |
Median and quantiles of specific columns |
median_hdi |
Median and Highest Probability Density Interval of specific columns |
median_hdci |
Median and continuous Highest Probability Density Interval of specific columns |
tidyMCMC |
Median/mean and quantiles/hpd of all columns |
It is important to remember that when working with Bayesian models, everything is a distribution. Whilst point estimates (such as a mean) of the parameters can be calculated from these distributions, we start off with a large number of estimates per parameter.
rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
fert.rstanarm3 %>% summary()
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: YIELD ~ FERTILIZER
## algorithm: sampling
## sample: 2400 (posterior sample size)
## priors: see help('prior_summary')
## observations: 10
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 52.7 13.0 36.7 52.7 68.9
## FERTILIZER 0.8 0.1 0.7 0.8 0.9
## sigma 19.1 4.9 13.8 18.2 25.5
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 163.4 7.9 153.4 163.5 173.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.3 1.0 2425
## FERTILIZER 0.0 1.0 2408
## sigma 0.1 1.0 2058
## mean_PPD 0.2 1.0 2443
## log-posterior 0.0 1.0 2338
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
Conclusions:
tidyMCMC(fert.rstanarm3$stanfit,
estimate.method = "median", conf.int = TRUE,
conf.method = "HPDinterval", rhat = TRUE, ess = TRUE
)
Conclusions:
fert.draw <- fert.rstanarm3 %>% gather_draws(`(Intercept)`, FERTILIZER, sigma)
## OR via regex
fert.draw <- fert.rstanarm3 %>% gather_draws(`.Intercept.*|FERT.*|sigma`, regex = TRUE)
fert.draw
We can then summarise this
fert.draw %>% median_hdci()
fert.rstanarm3 %>%
gather_draws(`(Intercept)`, FERTILIZER, sigma) %>%
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
Conclusions:
fert.rstanarm3 %>% plot(plotfun = "mcmc_intervals")
This is purely a graphical depiction on the posteriors.
fert.rstanarm3 %>% tidy_draws()
fert.rstanarm3 %>% spread_draws(`(Intercept)`, FERTILIZER, sigma)
# OR via regex
fert.rstanarm3 %>% spread_draws(`.Intercept.*|FERT.*|sigma`, regex = TRUE)
fert.rstanarm3 %>%
posterior_samples() %>%
as_tibble()
fert.rstanarm3 %>%
bayes_R2() %>%
median_hdci()
It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivariate distribution with elliptical contours.
mcmcpvalue <- function(samp) {
## elementary version that creates an empirical p-value for the
## hypothesis that the columns of samp have mean zero versus a
## general multivariate distribution with elliptical contours.
## differences from the mean standardized by the observed
## variance-covariance factor
## Note, I put in the bit for single terms
if (length(dim(samp)) == 0) {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp), transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1]) / length(samp)
}
else {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp), transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1]) / nrow(samp)
}
}
mcmcpvalue(as.matrix(fert.rstanarm3)[, c("FERTILIZER")])
## [1] 0
brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
fert.brm3 %>% summary()
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: YIELD ~ FERTILIZER
## Data: fert (Number of observations: 10)
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup samples = 2400
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 52.25 8.81 35.25 69.48 1.00 2470 1877
## FERTILIZER 0.81 0.06 0.69 0.92 1.00 2387 2287
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 12.78 1.87 9.62 16.84 1.00 2501 2240
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
fert.brm3$fit %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE,
conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE
)
Conclusions:
Return the draws (samples) for each parameter in wide format
fert.brm3 %>% tidy_draws()
The above can be useful for exploring the full posteriors of all the parameters of performing specific calculations using the posteriors, summarising the parameters takes a few more steps.
fert.brm3 %>%
tidy_draws() %>%
gather_variables() %>%
median_hdci()
The gather_draws on the other hand, conveniently combines tidy_draws and gather_variables together in a single command. Furthermore, it returns all of the variables. The spread_draws function allows users to target specific variables (either by naming them in full or via regexp).
fert.brm3 %>% get_variables()
## [1] "b_Intercept" "b_FERTILIZER" "sigma" "prior_Intercept"
## [5] "prior_b" "prior_sigma" "lp__" "accept_stat__"
## [9] "stepsize__" "treedepth__" "n_leapfrog__" "divergent__"
## [13] "energy__"
fert.brm3 %>%
gather_draws(b_Intercept, b_FERTILIZER, sigma) %>%
median_hdci()
## OR via regex
fert.brm3 %>%
gather_draws(`b_.*|sigma`, regex = TRUE) %>%
median_hdci()
fert.brm3 %>%
gather_draws(b_Intercept, b_FERTILIZER, sigma) %>%
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
fert.brm3 %>%
gather_draws(b_Intercept, b_FERTILIZER, sigma) %>%
ggplot() +
stat_halfeye(aes(x = .value, y = .variable))
Conclusions:
Select just the parameters of interest.
fert.brm3 %>% spread_draws(b_Intercept, b_FERTILIZER, sigma)
# OR via regex
fert.brm3 %>% spread_draws(`b_.*|sigma`, regex = TRUE)
fert.brm3 %>% mcmc_plot(type = "intervals")
Just the main parameters
fert.brm3 %>%
posterior_samples() %>%
as_tibble()
fert.brm3 %>%
bayes_R2(summary = FALSE) %>%
median_hdci()
It is possible to obtain an empirical p-value (for those that cannot live without p-values). This is essentially compares (for any posterior) that the mean is zero compared to a multivariate distribution with elliptical contours.
mcmcpvalue <- function(samp) {
## elementary version that creates an empirical p-value for the
## hypothesis that the columns of samp have mean zero versus a
## general multivariate distribution with elliptical contours.
## differences from the mean standardized by the observed
## variance-covariance factor
## Note, I put in the bit for single terms
if (length(dim(samp)) == 0) {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp), transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1]) / length(samp)
}
else {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp), transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1]) / nrow(samp)
}
}
mcmcpvalue(as.matrix(fert.brm3)[, c("b_FERTILIZER")])
## [1] 0
There are a large number of candidate routines for performing prediction. We will go through many of these. It is worth noting that prediction is technically the act of estimating what we expect to get if we were to collect a single new observation from a particular population (e.g. a specific level of fertilizer concentration). Often this is not what we want. Often we want the fitted values - estimates of what we expect to get if we were to collect multiple new observations and average them.
So while fitted values represent the expected underlying processes occurring in the system, predicted values represent our expectations from sampling from such processes.
| Package | Function | Description | Summarise with |
|---|---|---|---|
rstantools |
posterior_predict |
Draw from the posterior of a prediction (includes sigma) - predicts single observations | tidyMCMC() |
rstantools |
posterior_linpred |
Draw from the posterior of the fitted values (on the link scale) - predicts average observations | tidyMCMC() |
rstantools |
posterior_epred |
Draw from the posterior of the fitted values (on the response scale) - predicts average observations | tidyMCMC() |
tidybayes |
predicted_draws |
Extract the posterior of prediction values | median_hdci() |
tidybayes |
fitted_draws |
Extract the posterior of fitted values | median_hdci() |
tidybayes |
add_predicted_draws |
Adds draws from the posterior of predictions to a data frame (of prediction data) | median_hdci() |
tidybayes |
add_fitted_draws |
Adds draws from the posterior of fitted values to a data frame (of prediction data) | median_hdci() |
emmeans |
emmeans |
Estimated marginal means from which posteriors can be drawn (via tidy_draws |
median_hdci() |
For simple models prediction is essentially taking the model formula complete with parameter (coefficient) estimates and solving for new values of the predictor. To explore this, we will use the fitted model to predict Yield for a Fertilizer concentration of 110.
We will therefore start by establishing this prediction domain as a data frame to use across all of the prediction routines.
## establish a data set that defines the new data to predict against
newdata <- data.frame(FERTILIZER = 110)
fert.rstanarm3 %>%
fitted_draws(newdata = newdata) %>%
median_hdci()
fert.rstanarm3 %>% emmeans(~FERTILIZER, at = newdata)
## FERTILIZER emmean lower.HPD upper.HPD
## 110 141 129 152
##
## Point estimate displayed: median
## HPD interval probability: 0.95
fert.rstanarm3 %>%
posterior_predict(newdata = newdata) %>%
tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
fert.rstanarm3 %>%
posterior_epred(newdata = newdata) %>%
tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
fert.rstanarm3 %>%
posterior_linpred(newdata = newdata) %>%
tidyMCMC(estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval")
newdata %>%
add_predicted_draws(fert.rstanarm3) %>%
median_hdci()
newdata %>%
add_fitted_draws(fert.rstanarm3) %>%
median_hdci()
## Establish an appropriate model matrix
Xmat <- model.matrix(~FERTILIZER, data = newdata)
### get the posterior draws for the linear predictor
coefs <- fert.rstanarm3 %>%
posterior_samples(pars = c("(Intercept)", "FERTILIZER")) %>%
as.matrix()
fit <- coefs %*% t(Xmat)
fit %>% median_hdci()
fert.brm3 %>%
fitted_draws(newdata = newdata) %>%
median_hdci()
## OR
newdata %>%
add_fitted_draws(fert.brm3) %>%
median_hdci()
This is the option we will focus on for most of the worksheets
fert.brm3 %>%
emmeans(~FERTILIZER, at = newdata)
## FERTILIZER emmean lower.HPD upper.HPD
## 110 141 134 150
##
## Point estimate displayed: median
## HPD interval probability: 0.95
## OR
fert.brm3 %>%
emmeans(~FERTILIZER, at = newdata) %>%
tidy_draws() %>%
median_hdci()
fert.brm3 %>%
posterior_predict(newdata = newdata) %>%
as.mcmc() %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"
)
fert.brm3 %>%
posterior_epred(newdata = newdata) %>%
as.mcmc() %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"
)
fert.brm3 %>%
posterior_linpred(newdata = newdata) %>%
as.mcmc() %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"
)
newdata %>%
add_predicted_draws(fert.brm3) %>%
median_hdci()
newdata %>%
add_fitted_draws(fert.brm3) %>%
median_hdci()
## Establish an appropriate model matrix
Xmat <- model.matrix(~FERTILIZER, data = newdata)
### get the posterior draws for the linear predictor
coefs <- fert.brm3 %>%
posterior_samples(pars = c("b_Intercept", "b_FERTILIZER")) %>%
as.matrix()
fit <- coefs %*% t(Xmat)
fit %>% median_hdci()
Since we have the entire posterior, we are able to make probability statements. We simply count up the number of MCMC sample draws that satisfy a condition (e.g represent a slope greater than 0) and then divide by the total number of MCMC samples.
For this exercise, we will explore the following:
fert.rstanarm3 %>% hypothesis("FERTILIZER>0")
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (FERTILIZER) > 0 0.8 0.09 0.66 0.94 Inf 1
## Star
## 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Conclusions:
Alternatively…
fert.rstanarm3 %>%
tidy_draws() %>%
summarise(P = sum(FERTILIZER > 0) / n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 1
Conclusions:
newdata <- list(FERTILIZER = c(200, 100))
fert.rstanarm3 %>%
emmeans(~FERTILIZER, at = newdata) %>%
pairs()
## contrast estimate lower.HPD upper.HPD
## 200 - 100 80.5 63.5 97.4
##
## Point estimate displayed: median
## HPD interval probability: 0.95
Conclusions:
If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.
fert.mcmc <- fert.rstanarm3 %>%
emmeans(~FERTILIZER, at = newdata) %>%
tidy_draws() %>%
rename_with(~ str_replace(., "FERTILIZER ", "p")) %>%
mutate(
Eff = p200 - p100,
PEff = 100 * Eff / p100
)
fert.mcmc %>% head()
Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).
fert.mcmc %>% tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"
)
Alternatively, we could use median_hdci
fert.mcmc %>% median_hdci(PEff)
Conclusions:
To get the probability that the effect is greater than a 50% increase.
fert.mcmc %>% summarise(P = sum(PEff > 50) / n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 0.909
Conclusions:
Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).
fert.mcmc %>% hypothesis("PEff>50")
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0 10.7 8.35 -2.74 24.37 10.01 0.91
## Star
## 1
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.brm3 %>% hypothesis("FERTILIZER > 0")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (FERTILIZER) > 0 0.81 0.06 0.72 0.91 Inf 1
## Star
## 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.brm3 %>%
hypothesis("FERTILIZER > 0") %>%
plot(ignore_prior = TRUE)
fert.brm3 %>%
hypothesis("FERTILIZER > 0") %>%
`[[`("samples") %>%
median_hdci()
Conclusions:
Alternatively…
fert.brm3 %>%
gather_draws(b_FERTILIZER) %>%
summarise(P = sum(.value > 0) / n())
## # A tibble: 1 x 2
## .variable P
## * <chr> <dbl>
## 1 b_FERTILIZER 1
# OR
fert.brm3 %>%
tidy_draws() %>%
summarise(P = sum(b_FERTILIZER > 0) / n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 1
Conclusions:
fert.brm3 %>% hypothesis("FERTILIZER>0.9")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (FERTILIZER)-(0.9) > 0 -0.09 0.06 -0.18 0.01 0.06
## Post.Prob Star
## 1 0.06
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.brm3 %>%
hypothesis("FERTILIZER>0.9") %>%
plot(ignore_prior = TRUE)
# This returns a list
fert.brm3 %>%
hypothesis("FERTILIZER>0.9") %>%
plot(ignore_prior = TRUE) %>%
`[[`(1) +
geom_vline(xintercept = 0, linetype = "dashed")
# OR
fert.brm3 %>%
tidy_draws() %>%
summarise(P = sum(b_FERTILIZER > 0.9) / n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 0.0592
fert.brm3 %>%
tidy_draws() %>%
ggplot(aes(x = b_FERTILIZER)) +
geom_density(fill = "orange") +
geom_vline(xintercept = 0.9, linetype = "dashed")
newdata <- list(FERTILIZER = c(200, 100))
fert.brm3 %>%
emmeans(~FERTILIZER, at = newdata) %>%
pairs()
## contrast estimate lower.HPD upper.HPD
## 200 - 100 80.9 70.5 93.1
##
## Point estimate displayed: median
## HPD interval probability: 0.95
Conclusions:
If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.
It is in combination with emmeans() that tidy_draws() really comes into its own.
fert.mcmc <- fert.brm3 %>%
emmeans(~FERTILIZER, at = newdata) %>%
tidy_draws() %>%
rename_with(~ str_replace(., "FERTILIZER ", "p")) %>%
mutate(
Eff = p200 - p100,
PEff = 100 * Eff / p100
)
fert.mcmc %>% head()
Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).
fert.mcmc %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"
)
Alternatively, we could use median_hdci() to focus on a specific column.
fert.mcmc %>% median_hdci(PEff)
fert.mcmc %>%
ggplot() +
geom_density(aes(x = PEff))
Conclusions:
To get the probability that the effect is greater than a 50% increase.
fert.mcmc %>% summarise(P = sum(PEff > 50) / n())
## # A tibble: 1 x 1
## P
## <dbl>
## 1 0.980
Conclusions:
Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).
fert.mcmc %>% hypothesis("PEff>50")
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (PEff)-(50) > 0 10.91 5.61 1.87 20.4 50.06 0.98
## Star
## 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
fert.list <- with(fert, list(FERTILIZER = seq(min(FERTILIZER), max(FERTILIZER), len = 100)))
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) %>% as.data.frame()
head(newdata)
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
geom_point(data = fert, aes(y = YIELD)) +
geom_line() +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
scale_y_continuous("YIELD") +
scale_x_continuous("FERTILIZER") +
theme_classic()
## spaghetti plot
newdata <- emmeans(fert.rstanarm3, ~FERTILIZER, at = fert.list) %>%
gather_emmeans_draws()
newdata %>% head()
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
geom_line(aes(group = .draw), alpha = 0.01) +
geom_point(data = fert, aes(y = YIELD))
fert.grid <- with(fert, list(FERTILIZER = modelr::seq_range(FERTILIZER, n = 100)))
newdata <- fert.brm3 %>%
emmeans(~FERTILIZER, at = fert.grid) %>%
as.data.frame()
head(newdata)
ggplot(newdata, aes(y = emmean, x = FERTILIZER)) +
geom_point(data = fert, aes(y = YIELD)) +
geom_line() +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
scale_y_continuous("YIELD") +
scale_x_continuous("FERTILIZER") +
theme_classic()
## spaghetti plot
newdata <- emmeans(fert.brm3, ~FERTILIZER, at = fert.grid) %>%
gather_emmeans_draws()
newdata %>% head()
ggplot(newdata, aes(y = .value, x = FERTILIZER)) +
geom_line(aes(group = .draw), color = "blue", alpha = 0.01) +
geom_point(data = fert, aes(y = YIELD))